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1 Introduction 



Resonance is a common thread which runs through almost every branch of 
physics, without resonance we wouldn't have radio, television, music, etc. 
Resonance causes an object to oscillate, sometimes the oscillation is easy to 
see (vibration in a guitar string), but sometimes this is impossible without 
measuring instruments (electrons in an electrical circuit). A well-known ex- 
ample with Tacoma Narrows Bridge (at the time it opened for traffic in 1940, 
it was the third longest suspension bridge in the world) shows how disastrous 
resonances can be: on the morning of November 7, 1940, the four month old 
Tacoma Narrows Bridge began to oscillate dangerously up and down, tore 
itself apart and collapsed. Though designed for winds of 120 mph, a wind 
of only 42 mph caused it to collapse. The experts did agree that somehow 
the wind caused the bridge to resonate, and nowadays, wind tunnel testing 
of bridge designs is mandatory. 

Another famous example are the experiments of Tesla who studied in 1898 
experimentally vibrations of an iron column and noticed that at certain fre- 
quencies specific pieces of equipment in the room would start to jiggle. Play- 
ing with the frequency he was able to move the jiggle to another part of the 
room. Completely fascinated with these findings, he forgot that the column 
ran downward into the foundation of the building, and the vibrations were 
being transmitted all over Manhattan. The experiments had started sort of a 
small earthquake in his neighborhood with smashed windows, swayed build- 
ings, and panicky people in the streets. For Tesla, the first hint of trouble 
came when the walls and fioor began to heave . He stopped the experiment 
as soon as he saw police rushing through the door. 

The difference between resonances in a human made system and in some 
natural phenomena is very simple. We can change the form of a bridge and 
stop the experiment by switching off electricity but we can not change the 
direction of the wind, the form of the Earth atmosphere or the sizes of an 
ocean. What we can try to do is to predict drastic behavior of a real physi- 
cal system by computing its resonances. While linear resonances in different 
physical systems are comparatively well studied, to compute characteristics 
of nonlinear resonances and to predict their properties is quite a nontrivial 
problem, even in the one- dimensional case. Thus, the notorious Fermi-Pasta- 
Ulam numerical experiments with a nonlinear ID-string (carried out more 
then 50 years ago) are still not fully understood [2j. On the other hand, 
nonlinear wave resonances in continuous 2D-media like ocean, space, atmo- 
sphere, plasma, etc. are well studied in the frame of wave turbulence theory 
jl] and provide a sound basis for qualitative and sometimes also quantita- 
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tive analysis of corresponding physical systems. The notion of nonlinear 
wave interactions is crucial in the wave turbulence theory P]. Excluding 
resonances allows to describe a nonlinear wave system statistically, by wave 
kinetic equations and power-law energy spectra of turbulence [5J , and to ob- 
serve this behavior in numerical experiments [6]. Direct computations with 
Euler equations (modified for gravity water waves, [7]) show that the exis- 
tence of resonances in a wave system yield some additional effects which are 
not covered by the statistical description. The role of resonances in the evolu- 
tion of water wave turbulent systems has been studied profoundly by a great 
number of researchers. One of the most important conclusions (for gravity 
water waves) made recently in [8] is the following: " The four-wave resonant 
interactions control the evolution of the spectrum at every instant of time, 
whereas non-resonant interactions do not make any significant contribution 
even in a short-term evolution." 

The behavior of a resonant wave system can be briefly described [9] as follows: 
1) not all waves take part in resonant interactions, 2) resonantly interacting 
waves form a few independent small wave clusters, such that there is no 
energy flow between these clusters, 3) including some small but non-zero 
resonance width into consideration does not destroy the clusters. A model 
of laminated wave turbulence [10] allows to describe statistical and resonant 
regimes simultaneously while methods to compute resonances numerically 
are presented in (idea) and in (implementation). Our main purpose 
here is to study the possibilities of a symbolic implementation of these general 
algorithms using the computer algebra system Mathematica. 

The implemented software can be executed with local installations of Math- 
ematica and the corresponding method libraries; however, we have also de- 
veloped a Web interface that allows to run the methods from any computer 
in the Internet via a conventional Web browser. The implementation strat- 
egy is simple and based on generally available technologies; it can serve as a 
blueprint for other mathematical software with similar features. 

We take as our principal example the barotropic vorticity equation in a 
rectangular domain with zero boundary conditions which describes oceanic 
planetary waves, and show how : (a) to compute interaction coefficients of 
corresponding dynamical systems, (b) to solve resonant conditions, (c) to 
construct the topological structure of the solution set, and (d) to use the 
software via a Web interface over the Internet. A short discussion concludes 
the paper. 
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2 Mathematical Background 



Wave turbulence takes place in physical systems with nonlinear dispersive 
waves thatare described by evolutionary dispersive NPDEs. The role of the 
evolutionary dispersive NPDEs in the theoretical physics is so important 
that the notion of dispersion is used for a physical classification of PDEs into 
dispersive and non-dispersive. The well-known mathematical classification of 
PDEs into elliptic, parabolic and hyperbolic equations is based on the form 
of equations and can be applied to the second order PDEs on an arbitrary 
number of variables. On the other hand, the physical classification is based 
on the form of solutions and can be applied to PDEs of arbitrary order and 
arbitrary number of variables. In order to construct the physical classification 
of PDEs, two preliminary steps are to be made: 1) to divide all variables into 
two groups - time- and space-like variables {t and x correspondingly); and 2) 
to check that the linear part of the PDE under consideration has a wave-like 
solution in the form of Fourier harmonic 

ip{x, t) = A expi[kx — ut] 

with amplitude A, wave-number k and wave frequency u. The direct substi- 
tution of this solution into the linear PDE shows then that u is an explicit 
function on k, for instance: 

A + i'x + i^xxx = uj{k) = k- 5P. 

If u! as a function on k is real- valued and such that d'^uj/dk'^ ^ 0, it is 
called a dispersion function and the corresponding PDE is called evolutionary 
dispersive PDE. If the dimension of the space variable x is more that 

), is called the wave- vector and the dispersion 
function u = uj{k) depends on the coordinates of the wave-vector. This 
classification is not complementary to a standard mathematical one. For 
instance, though hyperbolic PDEs normally do not have dispersive wave 
solutions, the hyperbolic equation iptt — Ci'^^'xx — = has them. 

In the huge amount of application areas of NPDEs (classical and quantum 
physics, chemistry, medicine, sociology, etc.) a nonlinear term of the cor- 
responding NPDE can be regarded as small. This is symbolically written 
as 

L{^) = -eN{i,) (1) 

where L and are linear and nonlinear parts of the equation correspondingly 
and £ is a small parameter defined explicitly by the physical problem setting. 
It can be shown that in this case the solution ^ of ([T]) can be constructed 
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as a combination of the Fourier harmonics with amphtudes A depending on 
the time variable and possessing two properties formulated here for the case 
of quadratic nonlinearity: 



• PI The amplitudes of the Fourier harmonics satisfy the following sys- 
tem of nonlinear ordinary differential equations (ODEs) written for 
simplicity in the real form 



with coefficients a, being functions on wave-numbers; 

• P2 The dispersion function and wave-numbers satisfy the resonance 
conditions 



The transition form ([I]) to ([2]) can be performed by some standard methods 
(for instance, multi-scale method [1^) which also yields the explicit form of 
resonance conditions. 

Keeping in the mind our main problem - to find a solution of ([T]) - one 
has to take care of the initial and boundary conditions. This is done in the 
following way: the case of periodic or zero boundary conditions yields integer 
wave numbers, otherwise they are real. Correspondingly, one has to find all 
integer (or real) solutions of ([3]), substitute corresponding wave-numbers into 
the coefficients and then look for the solutions of ([2]) with given initial 
conditions. 

One can see immediately a big problem which appears as soon as one has to 
solve a NPDE with periodical or zero boundary conditions. Indeed, disper- 
sion functions take different forms, for instance, 

uj'^ = k^, up = + ak, uj'^ = k, uj = a/k, oj = m / n{n + 1) ■■■,etc. 

with k = (m, n), k = y/mP + n^ and a being a constant. This means that ([3]) 
corresponds to a system of Diophantine equations of many variables, nor- 
mally 6 to 9, with cumulative degrees 10 to 16. Those have to be solved 
usually for the integers of the order ~ 10^, which means that computa- 
tions has to be performed with integers of order 10^^ and more. Original 



A, 

A2 

A, 



axA^A'i 
012A1A2, 



(2) 




(3) 
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algorithms to solve these systems of equations have been developed based on 
some profound results of number theory [12] and implemented numerically 

m- 

Further on, an evolutionary dispersive NPDE with periodic or zero boundary 
conditions is called 3-term mesoscopic system if it has a solution of the form 

s s 

i' = '^i'i{^,t) = Ai exp iffcjfi - uji\ 

1=1 i=l 

and there exists at least one triple {Ai, Ai^, ^12} ^ {^i} such that PI and 
P2 keep true with some nonzero coefficients ai, a, 7^ = 1, 2, 3. 



3 Equations for Wave Amplitudes 
3.1 Method Description 

The barotropic vorticity equation describing ocean planetary waves has the 
form 

with boundary conditions 

■0 = for x = 0,L^; y = 0,Ly. 

Here /? is a constant called Rossby number, e is a small parameter and 
the Jacobean has the standard form 

J da db da db 

' dx dy dy dx 

First we give a basic introduction on how a PDE can be turned into a system 
of ODEs by a multi-scale method. Using operator notation, our problem (jl]) 
is viewed as a perturbed version of the linear PDE L{iIj) = 0. We pick a 
solution of this equation, say tpQ, which is a superposition of several waves 
ipj, i.e. ipo = ^jVj! each being a solution itself. To construct a solution 
of the original problem we make the amplitudes time-dependent. As the size 
of the nonlinearity in ([1]) is just of order e, the amplitudes will vary only on 
time-scales 1/e times slower than the waves. Hence we define an additional 
time- variable ti := called "slow time" to handle this time scale. So we 
look for approximate solutions of ([T]) that have the following form 

s 

tpo{t,ti,x) = ^Aj{ti)Lpj{x,t) 
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which for e = is an exact solution. The exact solution of the equation 
is written as power series in e around tp^, i.e. tp = YlT=o'^k£^- ■'■^ '^^^ 
computation it is truncated up to maximal order m which in our case is 
m = 1, i.e. 

ij(t,ti,x) = iJo(t,ti,x) + iJi(t,ti,x)e. 

Plugging ip[t, ti,x) one has to keep in mind that, since ti = et, we now have 
^ = ^ + £-^ due to the chain rule. Equations are formed by comparing the 
coefficients of e''. For k = this gives back the linear equation, but we keep 
the equation for = 1. In particular, for (jlj) we arrive at 

dt dx 

In order to ([2]), we have to get rid of all other variables. This is done by 
integrating against the yj/s, i.e. {.,ipj)L2(^Q^, and averaging over (fast) time, 

i.e. limT^oo Jq ■ dt. 



3.2 The Implementation 

This method was implemented in Mathematica with order m = 1 in mind 
only. So it won't be immediately applicable to higher orders without some 
(minor) adjustments. The ODEs are constructed done by the function 

ODESystemELC?/') , N(?/^) , 4), 

{xi,..,x„}, t, domain, jacobian, m, s, A, linwav, 
{Ai , . . , Ap} , paramvalues] . 

Basically this function takes the problem together with the solution of the 
linear equation as input and computes the list of ODEs for the amplitudes 
as output. Its arguments are in more detail: 

• L(ip) , Y^iip): Linear and nonlinear part of equation ([1]), each applied 
to a symbolic function parameter. Derivatives have to be specified with 
Dt instead of D and the nonlinear part has to be a polynomial in the 
derivatives of the function. 

• ip: symbol used for function in Liip) , U{ip) 
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• {xi, . . . ,x„}, t: list of symbols used for space- variables, and symbol 
for time- variable 

• domain: The domain on which the equation is considered has to be 
specified in the form {{xi ,ininxi ,maxxi} , ...» {x„,ininx„,maxx„}}, 
where the bounds on Xj may depend on xi , . . . , Xj_i only. 

• jacobian: For integration the (determinant of the) Jacobian must also 
to be passed to the function. This is needed in case the physical domain 
does not coincide with the domain of the variables above, it can be set 
to 1 otherwise. 

• m, s: maximal power of e and number of waves considered 

• A: symbol used for amplitudes 

• linwav: General wave of the linear equation is assumed to have sep- 
arated variables, i.e. (p{x,t) = Bi{xi)-...-Bn{xn) exp{i9{xi, ...,Xn,t)), 
and has to be given in the form 

{Bi(xi), B„(x„), 6'(xi , . . . ,x„,t)}. 

• {Ai, . . . ,\p}: list of symbols of parameters the functions in linwav 
depend on 

• paramvalues: For each of the s waves explicit values of the parameters 
{Ai,...,Ap} have to be passed as a list of s vectors of parameter values. 

ODESystem [linearpart_ , nonlinearpart_ , f un_Symbol , vars_List , 
t_Symbol , domain_List , j acobian_ , ord_Integer , num_Integer , 
A_Symbol , linwav_List ,parains_List ,paramvalues_List] := 
Module [{B,theta, eq,k} , 

eq = PerturbationEqns [linearpart ,nonlinearpart , 

fun,vars,t ,ord] ; 
eq = PlugInGenericWaveTuple[eq,fun,vars,t,A,B,theta,nuin] 

/. fun[l]->(0&) ; 
eq = Table [Resonaiice2 [eq, linwav, varSjt ,params , A, Bjtheta, 
num, paramvalues ,k] , 
{k.num}] ; 

Map [Integrate [Simplify [# , And@@ (Function [B , B [ [2] ] <B [ [1] ] < 

B [ [3] ] ] /Odomain) ] * j acobian , 
SequenceQOdomain] & , 

eq,{2}] 

] 
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Internally this function is divided into three subroutines briefly described 
below. 

3.2.1 Perturbation Equations, General Form 

The flrst of the subroutines is 

PerturbationEqns [L (■?/') , 'Ritp) , ip , {xi,...,x„}, t, m] . 

As mentioned before we approximate the solution of our problem by a polyno- 
mial of degree m in e. This subroutine works for arbitrary m. In the flrst step 
we construct equations by coefficient comparison. Additional time-variables 
will be created automatically and labeled t [1] , . . . ,t [m] . The output is a 
list of m-|- 1 equations corresponding to the powers 5°, ...,5™". The implemen- 
tation is quite straightforward. First set if] = YlT=o'^k(t,ti, ...,tm,xi, ..^Xn)s^ 
in ([1]), where tk = e^t, i.e. ^ = ^ + Yl^=i^^'^- Then extract the coefficients 
of on both sides and assemble the equations. Finally replace e^t by 

tfc again. 

PerturbationEqns [linearpart_ ,nonlinearpart_ , f un_Symbol , 
vars_List , time_Symbol , ord_Integer] : = 
Module [{i , j , e , eq} , 

eq = ((linearpart == -e*nonlinearpart) 

/. {fun->Sum[e"i*fun[i] [time,Sequence@@Table [e"j* 

time,{j jOrd}] ,Sequence@@ 
DeleteCases [vars , time] ] , 

{i,0,ord}]}); 

eq = (eq /. ((Dt[#, __]->0)& /(§ Join [vars , {time , e}] )) ; 
eq = (Equal(§(§#)& /@ 

Transpose [Take [CoefficientList [#, e] ,l+ord]& /@ 
(ListOOeq)] ; 
eq /. Table [e~j*time->time [j] jOrd}-] 
] 

3.2.2 Perturbation Equations, Given Linear Mode 

In step two we set ipoit, ti,x) = X]j=i t) described above. This 

is done by the function 

PluglnGenericWaveTuple [eq, ip , {xi,...,x„}, t, A, B, 6, s] 
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where the first argument is the output of the previous step. The symbols B 
and 9 have to be passed for labehng the shape and phase functions respec- 
tively. The output consists of two parts. The first part of the list formulates 
the assumption L{(pj) = explicitly for each of the waves. This is not used in 
subsequent computations, but is provided as a way to check the assumption. 
The second part of the hst is the equation corresponding to the coefficients 
of e from the previous step, with ipo as above. As the task of this step is so 
short the implementation does not need further explanation. 

PluglnGenericWaveTuple [eq_List , f un_Symbol , vars_List , 

t_Symbol , A_Symbol , B_Symbol , theta.Sjmibol , num_Integer] : = 
Module [{i , j , waves , n=Length [DeleteCases [vars , t] ] } , 
waves = Table [A [j] [Slot [2]]* 

Product [B[i] [j] [Slot [i+2] ] ,-[i,n}]* 
Exp[I*theta[j] [SequenceOQTable [Slot [i+2] , 

{i,n}] ,Slot[l]]] , 

{j .nuin}] ; 

{Table [eq[ [1] ] /. fun[0] ->Function[Evaluate [waves [ [j]]]] , 
{ j , num}] , 
Expand /O 

(eq[[2]] /. fun[0] ->Function [Evaluate [Total [waves] ]] ) 

}] 



3.2.3 Time and Scale Averaging 

Step three is the most elaborate. Under the assumption that interchange of 
averaging over time and inner product is justified, an integrand 

1 '•^ 



h= lim — I ipoipk dt 

T— >oo 1 







is computed that when integrated over the domain yields 

J ^^^^f J {'^o^Vk)L^n)dt. 
n ° 

Resonance conditions posed on the phase functions are explicitly used by 

Resonance [eq, linwav, {xi,..,x„}, t, 
{Ai,..,Ap}, A, B, 6, s, cond, k] 
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which receives the output from the previous step in eq. Here cond specifies 
the resonance condition in terms of the 6j, which have to be entered as 
6* [j] [xi , . . ,x„ ,t] respectively. The last argument is the index of the wave 
ipk in the integral above. Alternatively Resonance2 uses explicit parameter 
settings paramvalues for the waves instead of cond. This has been necessary 
because the general Resonance does not give useable results (see Section [X^ 
for more details). The main work in this step is to find out which terms do 
not contribute to the result. We exploit the fact that oscillating terms vanish 
when averaged over time by simply omitting those summands of {ipo, ipk) L'^{n) 
that have a factor exp(«6') with some time-dependent phase 9. The code for 
Resonance is not shown here, but is quite similar to Resonance2. 

Resonance2 [eq_List , linwav_List , vars_List , t_Symbol , params_List , 
A_Symbol , B_Symbol , theta_Symbol , num_Integer , 
paramvalues_List , testwave_Integer] : = 
Module [-[e , i , j , n=Length [DeleteCases [vars , t] ] } , 
e = Expand [ (List@@Last [eq] )* 

Exp [-l*theta [testwave] [Sequence®® 

DeleteCases [vars, t] , 
t]]] ; 

e = e / . 

Table [ 

theta[j] -> 

(Evaluate [(linwav [ [n+1] ] /. 

(RuleO{§#& /@ 
Transpose [{params , paramvalues [ [ j ] ] }] 

) 

) /. Append [Table [ 

DeleteCases [vars ,t] [[i]] 

-> Slot[i] , 
{i,n}] , 
t -> Slot [n+1]] 

]& 

), 

{ j , num}] ; 
e = MapAt [ 

(Function [theta, If [FreeQ [theta,t] ,theta,0] 
] 

[Simplify [#]] 

)&, 
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Position [e , Exp [_] ] ] ; 
e = Equal @@ 

(e*Conjugate [A [testwave] ] [t [1] ] * 
Product [Conjugate [B [i] 

[testwave] 

[DeleteCases [vars , t] [ [i] ] ] 

]. 

{i,n}] 

) /. 

Flatten [ 

Table [B[i] [j] -> 
Function [ 

Evaluate [DeleteCases [vars , t] [ [i] ] ] , 
Evaluate [linwav [ [i] ] /. 

(Rule(§@#& /(§ 
Transpose [ 

{params .paramvalues [ [j] ] 
>] 

)]], 

{i,n},{j ,num}]] 

] 

The integration of h is done by Mathematica and can be quite time-consuming. 
So ODESystem simplifies tlie integrand first to make integration faster. Still 
the expressions involved can be quite complicated. This is the most time- 
consuming part during construction of the ODEs. 



3.3 Obstacles 

Mathematica sometimes does not seem to take care of special cases and con- 
sequently has problems with evaluating expressions depending on symbolic 
parameters. We give two simple examples to illustrate this issue: 

• Orthogonality of sine-functions. 
Indeed, it holds that 

Vm, n e N : / sin(ma;) sm{nx)dx — 7r5jn,n- 
Jo 

Computing this in Mathematica by 
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Integrate [Sin [in*x] Sin [n*x] , {x , , 27r} , 

Assumptions —>■ mGlntegers && nGlntegers] 

yields independently of m, n instead. 

• Computation of a limit. 

Mathematica evaluates an expression 

, sinfxvr) 
Vn G Z : lim = Trdn o 

and similar expressions in two different ways getting two different an- 
swers. On the one hand 

Limit [Sin [(m-n)7r] /(m-n) , m^n, 

Assumptions mGlntegers && nGlntegers] 

gives 0. On the other hand, however, when the condition m,n E 1j 
is not used for computing the result Mathematica yields the correct 
answer vr, as with 

Limit [Sin [(m-n) tt] /(m-n) , m^n] . 

Unfortunately these issues prevented us from obtaining a nice formula for 
the coefficients in symbolic form by Resonance. So we just compute results 
for explicit parameter settings using Resonance2. 



3.4 Results 

3.4.1 Atmospheric Planetary Waves 

For the validation of our program we consider the barotropic vorticity equa- 
tion on the sphere first. Here numerical values of the coefficients are 
available (Table 1, [IS]). The equation looks quite similar 

However in spherical coordinates (0 G [— f , f ], A G [0,27r]) the differential 
operators are different: 

1 d 
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The linear modes have in this case the following form [Hj 



Pr (sin(0)) exp(.(mA + rrt)) (5) 

where P™(/i) are the associated Legendre polynomials of degree n and order 
m < n, so again they depend on the two parameters m and n. Also resonance 
conditions on the parameters look different in this case. 

Now we compute the coefficient as in ([2]). In [16] we find the following 
equation for the amplitude ^3 

OA 

nsins + l)^(^i) = 2z^K(n2 + 1) - n^in^ + l))Ai{t^)A2{h) 

so as = 2iZ . Parameter settings and corresponding numer- 
ical values for Z were taken from the table below (see [IS]). For this equa- 
tion and s = 3 results produced by our program have the form ciAs^s = 
C2A1A2A3, so as = C2/C1. 

Testing all resonant triads from the Table 1 from ^6], we see that the coeffi- 
cients differ merely by a constant factor of ±-\/8 which is due to the different 
scaling of the Legendre polynomials. In our computation they were nor- 
malized s.t. J^-^^ P^{fxYdfx = 1. With three triads, however, results were 
completely different. Interestingly this were exactly those triads for which 
no ipo appears in the table. 

Furthermore, for the other coefficients in ([2]) our program computes ai = 
a2 = in all tested parameter settings. This fact can be easily understood 
in the following way. We checked only resonance conditions but not the 
conditions for the interaction coefficients to be non-zero which are elaborated 
enough: 

nii < rii, rii 7^ rij Vi = 1, 2, 3, \ni — < n^, < rii + n2, 

and 

ni + n2 + ris is odd. 
Randomly taken parameter setting does not satisfy these conditions. 



15 



3.4.2 Ocean Planetary Waves 



Returning to the original example on the domain [0, L^:] x [0,Ly], we find 
explicit formulae for the coefficients in [I^. According to Section [X^ we can 
only verify special instances and not general formulae. 

Linear modes have now the form il5l 



sin(7r-— ) sin(7r— ^) exp(i(^x + cut)) (6) 

Ijr^ Ijy JjUJ 



with m, 77, G N and uj 



/3 



2-J(fr)^+{Tr)^ 



Parameter settings solving the resonance conditions were computed as in 
section HI Unfortunately results do not match and we have no explanation 
for that. In particular the condition ^ + + ^ = stated in [15] does not 
hold for the results of our program since we got ai = ^2 = in all tested 
parameter settings, just as in the spherical case. 

For example, if we try the triad {{2, 4}, {4, 2}, {1,2}} where = Ly = 1 
our program computes as = ^^j^vr (sin(3v^7r) — i{l + cos(3v^7r))) , whereas 
the general formula yields as = ^9+7\/5 ^ sin(3v^7r). However, if we use a triad 
with g = 1, e.g. {{24, 18} , {9 , 12} , {8 , 6}}, both agree on ai = a2 = as = 0. 



4 Resonance Conditions 

The main equation to solve is 



1 1 1 




for all possible mj, G Z with the scales and Ly (also G Z ) and then 
to check the condition rii ±n2 = n^. In the following argumentation it will 
be seen that and Ly can be assumed to be free of common factors. Below 
we refer to and Ly as to the scale coefficients. 

The first step of the algorithm implemented in Mathematica is to rewrite 
the equation to , ^ H — , ^ = , ^ and transform it in the 

yrni^+ni^ Vm2^+rr2^ Vma^+ris^ 

following way: we factorize the result of each rfij^ + and obtain with 
Pi - . . .■ Pr being the factors of mf + n1 and ai ■ . . . ■ a^ their respective powers: 

m^,+nj = pT-pT-...-pT- 
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We will now define a weight 7^ of the wave-vector {mi,ni) as the product 
of the pj's to the quotient of their respective aj and 2. The weight qi will be 
the name of the product of the p/s which have an odd exponent: 

Our equation then can be re-written as 

1 1 _ 1 

liy/qi l2y/q2 73v^ 

and one easily sees that the only way for the equation to possibly hold is 
Qi = ^2 = (l3 = Q (see [12] for details). Further we call q an index of the 
corresponding wave-vectors. The set of all wave-vectors with the same index 
is called a class of index q and is denoted as Clg. Obviously, the solutions 
of the resonance conditions are to be searched for with separate classes only. 

At this point one can also see that only such scales, and Ly, without 
common factors are reasonable. If they had a common factor, it would cancel 
out in the equation. 



4.1 Method Description 

The following five steps are the main steps of the algorithm: 

• Step 1: Compute the list of all possible indexes q. 

To compute the list of all indexes q, we use the fact that they have to 
be square-free and each factor of q has to be different from 3 mod 4 
(Lagrange theorem). There exist 57 possible possible indexes in our 
computational domains g < 300 : 

{1, 2, 5, 10, 13, 17, 26, 29, 34, 37, 41, 53, 58, 61, 65, 73, 74, 82, 85, 89, 
97, 101, 106, 109, 113, 122, 130, 137, 145, 146, 149, 157, 170, 173, 178, 
181, 185, 193, 194, 197, 202, 205, 218, 221, 226, 229, 233, 241, 257, 
265, 269, 274, 277, 281, 290, 293, 298} 

• Step 2: Solve the weight equation — + — = —. 

For solving the weight equation, we transform it into the equivalent 
form: 

73 = —— (7) 

7i + 72 
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The solution triples {71,72,73} can now be found by the two for-loops 
over 7i and 72 up to a certain maximum parameter and 73 is then being 
founded constructively with formula d?]). 

• Step 3: Compute all possible pairs {mi,ni) - if there are any - that 
satisfy mf + nf = •yf q. 

To compute our initial variables mi,ni, we use the Mathematica stan- 
dard function SumOfSquareRepresentation[d, x] which produces a 
list of all possible representations of an integer x as a sum of d squares, 
i.e. we can find all possible pairs (a, b) with d = 2 such that they 
satisfy + 6^ = x. Therefore, checking the condition mj + nj = •jf q 
is easy. 

• Step 4: Sort out the solutions {ml, ral, m2, ?t,2, m3, nS} that do not 
fulfill the condition nl ±n2 = n3. 

• Step 5: Check if by dividing the rrii by and the rii by Ly there are 
still exist some solutions. 

Last two steps are trivial. 
4.2 The Implementation 

Our implementation is quite straightforward and the main program is based 
on 4 auxiliary functions shown in the following subsections. 

4.2.1 List of Indexes 

The function constructqs[max] produces the list of all possible indexes q 
up to the parameter max. The first (obvious) g's sol = {1} is given and the 
function checks the conditions starting with n = 2. Every time n satisfies 
the conditions, it is appended to the list sol. If one condition fails, the next 
n = n + 1 is considered and so on until n reaches the parameter max. Then 
the list sol is returned: 

Clear [constructqs] ; 

constructqs [n_, sol_List, max_] ; n>max := sol (*6*) 

constructqs [n_?SquareFreeQ , sol_List, max_] 

:= constructqs [n+1 , Append [sol, n] , max] (*5*) 



18 



constructqs [n_?SquareFreeQ , sol_List, maxJ] ; 
MemberQ[Mod[PrimeFactorList [n] , 4], 3] 
:= constructqs [n+1 , sol, max] (*4*) 

constructqs [n_, sol_List, max_] ; ! SquareFreeQ [n] 
:= constructqs [n+1 , sol, max] (*3*) 

constructqs [1] := {1} (*2*) 

constructqs [max_] := constructqs [3, {l}, max] (*1*) 
4.2.2 Weight Equation 

The function find7s[7max] solves the weight equation in the following way. 
For a fixed 71 and 72 running between 1 and ■ymax, it is checked if 73 is 
an integer. If it is, the triple {71,72,73} is added to the list sol which is 
empty at the initial moment. Once 72 reaches ^ymax, it is set to 1 again 
and the search starts again with 71 = 71 + 1. This is done as long as both 
7i and 72 are lower than max. Finally the list sol is returned: 

f ind7S [7max;_, 7l_, j2_, sol_List] ; 

7I > 7max := (Clear [73] , sol) (*6*) 

f ind7S [7max;_ , 7l_, 72_, sol_List] ; (7I < 7max && 72>7max; && 
IntegerQ [73= (7I72) / (7I+72) ] ) 

:= f ind7S [7max; , 7I+I, 1, Append[sol, {7I, 72, 73}]] (*5*) 

f ind7s [7max_, 7I-, 72_, sol_List] ; 

(7I < 7max && 72>7max && 

! IntegerQ [73= (7172)7(71+72)] ) 

:= find7s[7max, 7I + 1, 1, sol] (*4*) 

f ind7S [7max_, 7l_, 72_, sol_List] ; 

(7I < 7max && 72 < 7max && IntegerQ [73= (7172)7(71+72) ] ) 

:= f ind7S [7max;, 7I, 72 + 1, Append [sol, {7I, 72, 73}]] (*3*) 

f ind7S [7max;_, 7l_, 72_, sol_List] ; 

(7I < 7max && 72 < 7max && ! IntegerQ [73= (7I72) 7(71+72) ] ) 
:= find7s[7max, 7I, 72 + 1, sol] (*2*) 

f ind7S [7max;_] := find7s[7max, 1, 1, {}]) (*1*) 

For flnd7s[7max] to be executable, the iteration depth of 2^^ is not sufficient 
and it was set to 00. 
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4.2.3 Linear Condition 



The third auxihary function makemns checks whether the hnear condition 
ni ± 77,2 = ^^3 is fulfilled and structures the solution set into a list of pairs 

{{mi, ni}, {m2, 712}, {7773, 773}} : 

Clear [makemns] ; 

makemns [ml_, nl_, m2_, n2_, m3_, n3_] := {} (*3*) 

maiemns [ml_, nl_, m2_, n2_, m3_, n3_] ; 
(nl + n2 == n3 || nl - n2 == n3) : = 

{{ml, nl}, {m2, n2}, {m3, n3}} (*2*) 

makemns [mnl_List , mn2_List, mn3_List] : = 

Cases [Flatten [Table [makemns [mnl [[1,1]] , mnl [[1,2]] , 
mn2[[j,l]], mn2[[j,2]], mn3[[k,l]], mn3[[k,2]]], 
{i, 1, Length [mnl] } , { j , 1, Length [mn2] } , 
{k, 1, Length [mn3] }] , 2], 
{{xl_,x2_}, {x3_,x4_}, {x5_,x6_}}] (*1*) 

The function makemns is called three times: 

In from 3 lists of arbitrarily many pairs {mi, ni}, a 3-dimensional 

array is made combining entries of the 3 lists with each other. Each entry 
calls the same program with the parameters of the current combination of 
{ml,nl,m2,n2,m3,n3}. 

In (*2*) and (*3*) it is decided whether the condition 77I ± n2 = n2> 
is fulfilled. If it is, a solution {{ml,nl},{m2,n2},{m3,n3}} is written in 
the array. The table is then flattened to the level 2 in order to have a 
list of solutions. In the end, all empty lists have to be sorted out, done 
by the function Ceises which keeps only those cases that have the shape 
{{xl_,x2_},{x3_,x4_},{x5_,x6_}}. 

4.2.4 Scale Coefficients 

Finally, the function respectL[soI, Lx, Ly] divides each component of the 
solution by the pair (L^,, Ly) and sorts out the result if any of the 6 compo- 
nents does not remain an integer: 

respectL [sol_Llst , Lx_, Ly_] : = 
Map [solution [#] k , 

Cases [Map [#/{Lx, Ly}&, 

Map[#[[l]]]&, sol], {2}], {{.Integer, .Integer}, 
{-Integer, .Integer}, {.Integer, .Integer}}]] 
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The function respectL[sol, Lx, Ly] gets as an input the hst of the form 
{solution[{{ml,nl},{m2,n2},{m3,n3}}],...} and returns the hst of the same 
form. 

4.3 Results 

All solutions in the computation domain m,n < 300 have been found in 
a few minutes. Notice that computations in the domain m,n < 20 by 
direct search, without introducing indexes q and classes Clq took about 
30 minutes. A direct search in the domain m,n < 30 has been interrupted 
after 2 hours, since no results were produced. 

The number of solutions depends drastically on the scales and Ly, 
some data are given below (for the domain m, n < 50 : ) 

{L^ = 1, Ly = 1) : 76 solutions; 

(La; = 3, Ly = 1) : 23 solutions; 

{Lx = 6, Ly = 16) : 2 solutions; 

{Lx = 5, Ly = 21) : 2 solutions; 

{Lx = II, Ly = 29) : no solutions (search up to 300, for both qmax and 
'jmax) . 

Interestingly enough, in all tried possibilities, only an odd q yield solutions. 

5 Structure of the Solution Set 
5.1 Method Description 

The graphical way to present 2D-wave resonances suggested in [9j for 3- wave 
interactions is to regard each 2D- vector k = {m, n) as a node (m, n) of inte- 
ger lattice in the spectral space and connect those nodes which construct one 
solution (triad, quartet, etc.). Having computed already all the solutions of 
([2D in Section 4, now we are interested in the structure of resonances in spec- 
tral space. To each node (m, n) we can prescribe an amplitude A{m, n, ti) 
whose time evolution can be computed from the dynamical equations ob- 
tained in Section 3. Thus, solution set of resonance conditions ([3]) can be 
thought of as a collection of triangles, some of them are isolated, some form 
small groups connected by one or two vertices. Corresponding dynamical 
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systems can be re-constructed from the structure of these groups. For in- 
stance, a single isolated triangle corresponding to a solution with wave vectors 
(mi, ni)(m2, n2){ms, n^) and wave amplitudes {{Al, A2, A3)} corresponds to 
the following dynamical system: 

A2 = 0:2^1^3 
A3 = a3AiA2 

with being functions of all rrii, rii (see Section 3). 

If that two triangles share one common vertex {(741,742,743), (743,744,745)}, 
the the corresponding dynamical system is 

74i = aiA2A^ 
A2 = a2AiAs 

A3 = ^3 J. 4 1.42 + a3^2A4A5 

A4 = a^A'sA^ 
A5 = a^A^Ai 

If two triangles have two vertices in common {(741, 742, 743), (742, 743, 744)}, 
then the dynamical system is quite different: 

Al = aiA2A3 

A2 = a2,iAiA3 + 0:2,2^3^4 

743 = Q;3^174i742 + (^3,2^2^4 

A4 = aiA2A3 = — ii 

Using the fourth equation, the formulae for .42 and .43 can be simplified to: 

744 = —Al =^ 744 = — 74i (5i 

ai ai 

A2 = 74i743 Qi2,i H ■ H 

V «! y «! 

A3 = 74i742 q;3,i H ■ H 

V ai J ai 

This means that qualitative dynamics of the 3-term mesoscopic system de- 
pends not on the geometrical structure of the solution set but on its topological 
structure. Constructing the topological structure of the solution set, we do 
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not consider concrete values of the solution but only the way how triangles 
are connected. In any finite spectral domain we can compute all independent 
wave clusters and write out corresponding dynamical systems thus obtain- 
ing complete information about energy transfer through the spectrum. Of 
course, quantitative properties of the dynamical systems depend on the spe- 
cific values of rrii^ni (for instance, values of interaction coefficients ctj, 
magnitudes of periods of the energy exchange among the waves belonging to 
one cluster, etc.) 

5.2 Implementation 

To construct the topological structure of a given solution set we need first to 
find all groups of connected triangles. This is done by the following proce- 
dure: 

FindConnectedGroups [triangles_List] := 
Block [{groups = {}, tr = triangles, newgroup}, 
While [Length [tr] > 0, 
{newgroup, tr} = 

FindConnectedTriangles [{First [tr] } , Rest [tr] ] ; 
groups = Append [groups , newgroup] ; 

]; 

groups 

]; 

FindConnectedTriangles [grp_List , triangles_List] : = 
Module [{points , newGrpMember , tr=triangles} , 
points=Flatten [Apply [List ,grp, 2] ,1] ; 

newGrpMember=Cases [tr , _ [ , #1 , ] ] &/ ©points ; 

(tr=DeleteCases [tr , _ [ , #1 , ] ] )&/@points ; 

newGrpMember=Union [JoinQOnewGrpMember] ; 
If [Length [newGrpMember] ==0 , 
{grp,tr}, 

newGrpMember=FindConnectedTriangles [newGrpMember , tr] ; 
{Join [grp , First [newGrpMember] ] , 
newGrpMember [ [2] ] } 

] 

]; 

The function FindConnectedGroups expects a list of triangles as input, and 
three different types for data structure can be used. The first type is just 
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a list of three pairs, where each pair contains the coordinates of a node, for 
example {{1,2}, {3, 4}, {5, 6}}. An alternative type is like the type before 
just with another head symbol instead of list, e.g. 
Triangle [{1 , 2} , {3 , 4} , {5 , 6}] . 

The function also works for vertex numbers instead of coordinates, e.g. 
Triangle [1, 2, 3]. In every case the function returns a partition of the 
input list where all elements of a list are connected and elements of different 
lists have no connection to each other. 

The function FindConnectedTriangles is an auxiliary function which has 
two parameters. The first list contains allconnected triangles. The second 
list contains all other triangles which are possibly connected to one of the 
triangles in the first list. The function FindConnectedTriangles returns a 
pair of lists: the first list contains all triangles which are connected to the 
selected triangles, the second list contains all remaining. 

The input list for FindConnectedTriangles is a list of 3-element lists. Before 
we can use the results produced in Section [Has an input we have to transform 
the data. This can be easily done by: 

Transf ormSolution [sol_List] := 

Flatten [Rest/Osol] / . solution [trs : { List}] ->trs ; 

Some remarks on the implementation. 

The function FindConnectedGroups selects a triangle, which is not yet in a 
group and calls the function FindConnectedTriangles. Since the returned 
first list always contains at least one triangle, the length of the list tr de- 
creases in every loop call, hence the FindConnectedGroups terminates. The 
question left is how to find all triangles connected with a certain triangle. 
This has been done in the following way. First we search for all triangles 
which share at least one node with this triangle. Then we restart the search 
with all triangles found. For efficiency reasons it is better to perform the 
search with all triangles we found in one step together. If in one step no fur- 
ther triangles are found then we are ready and return the list of connected 
triangles and the remaining list. In each step we remove all triangles we 
found from the list of triangles which are not declared as connected. This 
increases the speed because the search is faster if there are less elements to 
compare. More important, this prevent us to search in loops and find some 
triangles more than once. In general, search in a loop can be the reason for 
a termination problem but due to shrinking the list of triangles to search for 
in every step the termination can be guaranteed. 
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5.3 Results 

In the Figured] the geometrical structure of the solution set is shown, for the 
case mi,ni < 50 and = Ly = 1. 




Figure 1: The geometrical structure of the result in domain D = 50 

Below we show all the topological elements of this solution set. 

1. 21 groups contain only one triangle (obviously, they have isomorphic 
dynamical systems): 

{{3, 18}, {36, 6}, {2, 12}} {{4, 46}, {14, 44}, {23, 2}} 

{{6, 44}, {36, 26}, {13, 18}} {{6, 48}, {42, 24}, {3, 24}} 

{{8, 26}, {16, 22}, {13, 4}} {{9, 24}, {48, 18}, {16, 6}} 

{{14, 28}, {28, 14}, {7, 14}} {{18, 36}, {36, 18}, {9, 18}} 

{{22, 16}, {26, 8}, {11, 8}} {{22, 20}, {28, 10}, {11, 10}} 

{{22, 44}, {44, 22}, {11, 22}} {{22, 48}, {42, 32}, {21, 16}} 

{{24, 18}, {9, 12}, {8, 6}} {{26, 28}, {28, 26}, {19, 2}} 

{{28, 42}, {42, 28}, {21, 14}} {{28, 46}, {50, 20}, {7, 26}} 

{{36, 22}, {42, 4}, {11, 18}} {{36, 30}, {15, 18}, {10, 12}} 

{{38, 24}, {42, 16}, {21, 8}} {{44, 18}, {46, 12}, {23, 6}} 
{{48, 36}, {18, 24}, {16, 12}} 
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2. Further 9 groups contain also one triangle, but in each triangle two points 
coincide (again, they have isomorphic dynamical systems): 

{{8,2}, {8,2}, {1,4}} {{16,2}, {16,2}, {7,4}} 

{{16, 4}, {16, 4}, {2, 8}} {{24, 6}, {24, 6}, {3, 12}} 

{{32, 8}, {32, 8}, {4, 16}} {{34, 8}, {34, 8}, {7, 16}} 

{{46, 8}, {46, 8}, {17, 16}} {{48, 6}, {48, 6}, {21, 12}} 
{{48, 12}, {48, 12}, {6, 24}} 



10 20 



3. There exist 2 groups with two triangles each (by observation of the geo- 
metrical pictures it is easy to determine that both have isomorphic dynamical 
systems) : 

{ {{2, 24}, {18, 16}, {9, 8}}, {{4, 48}, {36, 32}, {18, 16}} } 
{ {{12, 26}, {26, 12}, {3, 14}}, {{26, 12}, {28, 6}, {13, 6}} } 





4. Two further groups consist of two triangles each, but the common point is 
contained twice in one triangle (the dynamical systems are isomorphic, but 
different from the two groups above): 

{ {{24, 22}, {32, 6}, {3, 16}}, {{32, 6}, {32, 6}, {11, 12}} } 
{ {{8, 38}, {32, 22}, {11, 16}}, {{38, 8}, {38, 8}, {11, 16}} } 
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5. As we can see by inspecting their geometrical structures, further 7 groups 
are not isomorphic to any group found above: 



{ 



{ 



{ 



{{6, 12}, {12, 6}, {3, 6}}, {{12, 24}, {24, 12}, {6, 12}}, 
{{24, 48}, {48, 24}, {12, 24}} 

{ {{2, 16}, {14, 8}, {1,8}}, {{4, 32}, {28, 16}, {2, 16}}, 
{{32, 4}, {32, 4}, {14, 8}} } 

{ {{2,4}, {4,2}, {1,2}}, {{4,8}, {8,4}, {2,4}}, 

{{8, 16}, {16, 8}, {4, 8}}, {{16, 32}, {32, 16}, {8, 16}} } 

{{4, 22}, {10, 20}, {11, 2}}, {{8, 44}, {20, 40}, {22, 4}}, 
{{10, 20}, {20, 10}, {5, 10}}, {{20, 40}, {40, 20}, {10, 20}} 

{{10, 40}, {26, 32}, {19, 8}}, {{26, 32}, {38, 16}, {13, 16}}, 
{{32, 26}, {40, 10}, {13, 16}}, {{40, 10}, {40, 10}, {5, 20}} 

{{4, 18}, {14, 12}, {7, 6}}, {{8, 36}, {28, 24}, {14, 12}}, 



{{12, 14}, {14, 12}, {9, 2}}, 
{{36, 42}, {42, 36}, {27, 6}}, 

{{2, 36}, {20, 30}, {17, 6}}, 
{{8, 12}, {12,8}, {6,4}}, 
{{16, 24}, {24, 16}, {12, 8}}, 



{{24, 28}, {28, 24}, {18, 4}}, 
{{42, 36}, {21, 18}, {4, 18}} 

{{4, 6}, {6, 4}, {3, 2}}, 
{{12, 18}, {18, 12}, {9, 6}}, 
{{18, 12}, {9, 6}, {4, 6}}, 



} 



{{20, 30}, {30, 20}, {15, 10}}, {{20, 30}, {34, 12}, {1, 18}}, 

{{24, 36}, {36, 24}, {18, 12}}, {{30, 20}, {36, 2}, {1, 18}}, 

{{32, 48}, {48, 32}, {24, 16}}, {{34, 12}, {36, 2}, {15, 10}}, 

{{36, 24}, {18, 12}, {8, 12}}, {{45, 30}, {34, 12}, {12, 18}} 
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5.4 Important Remark 



To compute all non-isomorphic sub-graphs algorithmically is a nontrivial 
problem. Indeed, all isomorphic graphs presented in previous section are 
described by similar dynamical systems, only magnitudes of interaction co- 
efficients ai vary. However, in the general case graph structure thus defined 
does not present the dynamical system unambiguously. Consider Figure [2] 
below where two objects are isomorphic as graphs. However, the first object 
represents 4 connected triads with dynamical system 

(/ll,A2,A3), (/li,A2,A5), (^1,^3,^4), (^2,^3,^6) (8) 

while the second - 3 connected triads with dynamical system 

(^1,^2,^5), (Ai,/l3,A4), (^2,^3,^6). (9) 




Figure 2: Example of isomorphic graphs and non-isomorphic dynamical sys- 
tems. The left graph corresponds to the dynamical system ([8]) and the graph 
on the right - to the dynamical system Q. To discern between these two 
cases we set a placeholder inside the triangle not representing a resonance. 



This problem has been solved in pLTj by introducing hyper-graphs of a special 
structure; the standard graph isomorphism algorithm used by Mathematica 
has been modified in order to suit hyper-graphs. 

6 A Web Interface to the Software 

The previous sections have presented implementations of various symbolic 
computation methods for the analysis of non-linear wave resonances. These 
implementations are written in the language of the computer algebra system 
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Mathematica which provides an appeahng graphical user interface (GUI) 
for executing computations and presenting the results. For instance, the 
pictures shown in Section 14.31 were produced by converting the computed 
hyper-graphs to Mathematica plot structures that can be displayed by the 
GUI of the system. 

However, to run these methods the user needs an installation of Mathematica 
on the local computer with the previously described methods installed in a 
local directory. These requirements make access to the software difficult and 
hamper its wide-spread usage. In order to overcome this problem, we have 
implemented a Web interface such that the software can be executed from 
any computer connected to the Internet via a Web browser without the need 
for a local installation of mathematical software. 

This implementation follows a general trend in computer science which turns 
away from stand alone software (that is installed on local computers and 
can be only executed on these computers via a graphical user interface) and 
proceeds towards service- oriented software [18] (that is installed on remove 
server computers and wraps each method into a service that can be invoked 
over the Internet via standardized Web interfaces). Various projects in com- 
puter mathematics have pursued middleware for mathematical web services, 
see for instance p^ l20 | [2T]. On the long term, it is thus envisioned that math- 
ematical methods generally become remote services that can be invoked by 
humans (or other software) without requiring local software installations. 

However, even without sophisticated middleware it is nowadays relatively 
simple to provide (for restricted application scenarios) web interfaces to 
mathematical software by generally available technologies. The web interface 
presented in the following sections is deliberately kept as simple as possible 
and makes only use of such technologies; thus it should be easy to take this 
solution as a blueprint for other mathematical software with similar features. 
In particular, the web interface is quite independent of Mathematica as the 
system underlying the implementation of the mathematical methods; the 
same strategy can be applied to other mathematical software systems such 
as Maple, MATLAB, etc. 

6.1 The Interface 

Figure [3] shows the web interface to some of the methods presented in the 
previous sections. Its functionality is as follows: 

Create Solution Set The user may enter a parameter D in the first (small) 
text field and then press the button "Create Solution Set" . This invokes 
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Hint: for plotting the t'Opobgy, you first need to 



a solution set 60111 a domain vake ct pro«de ■ 




Figure 3: Web interface to the implementation 

the method CreateSolutionSet which computes the set of all solutions 
whose values are smaller than or equal to D. This set is written into 
the second (large) text field in the form 

{Solution [xi ,2;i] , . . . , Solution [x„ } 

Plot Topology The user may enter into the second (large) text field a spe- 
cific solution set (or, as show above, compute one), and then press 
the button "Plot Topology". This first invokes the method Topology 
which computes the topological structure of the solution set as a list 
of hyper-graphs and then calls the method PlotTopology which com- 
putes a plot of each hyper-graph. The results are displayed in the right 
frame of the browser window. 

The web interface is available at the URL 

http : //www . rise . uni-linz . ac . at /pro j ects/ alisa 

(Button "Discrete Wave Turbulence") 

To run the computations, an account and a password are needed. 
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6.2 The Implementation 



The web interface is implemented in PHP, a scripting language for producing 
dynamic web pages [22j. PHP scripts can be embedded into conventional 
HTML pages within tags of form <php? . . . ?>; when a Web browser requests 
such a page, the Web server executes the scripts with the help of an embedded 
PHP engine, replaces the tags by the generated output, and returns the 
resulting HTML page to the browser. With the use of PHP, thus programs 
can be be implemented that run on a web server and deliver their results to 
a client computer which displays them in a web browser. The web interface 
to the discrete wave turbulence package is implemented in PHP as sketched 
in Figure H] and described below (the parenthesized numbers in the text refer 
to the corresponding numbers in the figure). 

Create Solution Set The browser frame input on the left side contains 
essentially the following HTML input form: 

<forni target="textarea" 

action="https : //apache2 . . . /CreateSolutionSet .php" 
method="post"> 
<input name=" domain" size="3"> 

<input type=" submit" value="Create Solution Set"> 
</f orm> 

This form consists of an input field domain to receive a domain value and 
a button to trigger the creation of the solution set. When the button is 
pressed, (1) a request is sent to the web server which carries the value of 
domain; this request asks the server to deliver the PHP-enhanced web page 
CreateSolutionSet .php into the target frame textarea which is displayed 
internally to input. 

The file CreateSolutionSet .php has essentially the content 

<?php 

$math="/. . ./math" ; 

$cwd="/. . . /DiscreteWaveTurbulence" ; 
$domain = $_POST[Momain'] ; 
$mcmd = 

"SetDirectory[\"" . $cwd . "\"] ; " . 
"Needs [\"DiscreteWaveTurbulence ' SolutionSet ' \"] ; " . 
"sol=DiscreteWaveTurbulence'SolutionSet 'CreateSolutionSet [" . 
$domain . "] ; " ; 
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Figure 4: Implementation of the web interface 
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$coinmaiid="$math -noprompt -run '" . $mcmd . 

"Print [StandardForm [sol] ] ; Quit □ ; ' " ; 
$result = shell_exec("$coinmaiid") ; 
echo 

"<textarea name=\"sol\" cols=\"60\" rows=\"20\">" . 
htmlspecial chars ($result) . 
"</textarea>" . 

. . . , 

?> 

After setting the paths $math of the Mathematica binary and $cwd of the 
directory where the DiscreteWaveTurbulence package is installed, the script 
sets the local variable Sdomain to the value of the input field domain. Then 
the Mathematica command Smcind is constructed in order to load the file 
SolutionSet .m and execute the command CreateSolutionSet to compute 
the solution set. Now the system command $coinmand is constructed to (2) 
invoke Mathematica which calls the previously constructed command and 
(3) prints its result to the standard output stream which is captured in the 
variable $result. Prom this, the script contstructs the HTML code of the 
result document which is (4) dehvered to the Web browser. 

Plot Topology The browser frame textarea contains essentially the fol- 
lowing HTML input form: 

<form target=" result" 

action="https://apache2. . . ./PlotTopology.php" 

method="post"> 
<textarea name="sol" cols="60" rows="20"> . . . </textarea> 
<input type="submit" value="Plot Topology"> 
</center> 
</f orm> 

This form consists of the textarea field sol to receive the solution set and a 
button to trigger the plotting of the topology of this set. When the button 
is pressed, (1) a request is sent to the web server which carries the value 
of sol; this request asks the server to deliver the PHP-enhanced web page 
PlotTopology.php into the target frame result on the right side of the 
browser. 

The file CreateSolutionSet .php has essentially the content 
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<?php 

$math="/. . ./math" ; 

$basedir ="/. . . /DiscreteWaveTurbulence" ; 

$baseurl ="http://apache2/. . ./DiscreteWaveTurbulence"; 

$sol = $_POST['sol'] ; 

... // create under $basedir a unique subdirectory $dir 
$mcmd = 

"SetDirectory [\"$basedir/$dir\"] ; " . 
"Needs [\"DiscreteWaveTurbulence 'Topology ' \"] ; " . 
"Needs[\"DiscreteWaveTurbulence'SolutionSet'\"] ; " . 
"top=DiscreteWaveTurbulence 'Topology' Topology [$sol] ; " . 
"plots=DiscreteWaveTurbulence' Topology 'PlotTopologyl [top] ; " ; 
$command="/usr/bin/Xvnc :20 & export DISPLAY=:20; " . 
"export MATHEMATICA_USERBASE=$basedir/.Mathematica; " . 
"$math -run '" . $mcmd . 

"Print [ExportList [plots, \"$image\"]] ; Quit [] ; '"; 
$result = shell_exec("$command | tail -n 1"); 
for ($i=0;$i<$result;$i++) 

echo "<img src=\"$baseurl/$dir/image-$i .png\"/>" ; 

?> 

For holding the images to be generated later, the script creates a unique 
directory $basedir/$dir which is served by the web server under the url 
$baseurl/$dir. The script extracts the solution set $sol from the request 
and sets up the Mathematica command to compute its topological structure 
and generate the plots from which ultimately the image files will be produced. 

For this purpose, however, Mathematica needs an Xll display server running; 
since a Web server has not access to an Xll server, we start the virtual Xll 
server Xvnc [23] as a replacement and set the environment variable DISPLAY 
to the display number on which the number listens; Mathematica will subse- 
quently send Xll requests to that display which will be handled by the virtual 
server. Likewise, Mathematica needs access to a .Mathematica configuration 
directory; the script sets the environment variable MATHEMATICA_USERBASE 
correspondingly. 

With these provisions, we can (2) invoke first the command to compute 
the plots and then the (self-defined) command ExportList to generate for 
every plot an image in the previously created directory. For this purpose the 
command uses (3) the Mathematica command EXPORT [/iie, plot, "PNG"] 
which converts plot to an image in PNG format and writes the image to file. 
ExportList returns the number of images generated which is (4) written 
to the standard output stream which in turn is captured in the variable 
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Sresult. From this information, the script generates an HTML document 
which contains a sequence of img elements referencing these images. After 
this document has been (5) returned to the chent browser, the browser (6) 
requests the referenced images with GET messages from the web server. 

6.3 Extensions 

As an alternative to the display of static images, the Web interface also 
provides an option "Applet Viewer" with somewhat more flexibility. If this 
option is selected, Mathematica is instructed to save all generated plots as 
files in the standard representation. The generated HTML document then 
embeds (rather than img elements) a sequence of applet elements that load 
instances of the "JavaView" applet [23]. These applets run in the Java Vir- 
tual Machine of the Web browser on the client computer, load the plot files 
from the web, and visualize them in the browser. Rather than just displaying 
static images, the viewer allows to perform certain manipulations and trans- 
formations of the plots such as scaling, rotating, etc. While this additional 
flexibility is not of particular importance for the presented methods, they 
may in the future become useful for others. 

To limit access to the software respectively to the computing power of the 
server computer, it may be protected by authentication mechanisms. For 
example, on the Apache Web server, it suffices to provide in the installation 
directory of the software a file .htaccess with content 

<Files "*.php"> 

SSLRequireSSL 

AuthName "your account" 

AuthType Basic 

Require valid-user 
</Files> 

With this configuration, the user is asked for the data of a valid account 
on the computer running the Web server; other authentication mechanisms 
based e.g. on password files may be provided in a similar fashion. 

7 Discussion 

Summing up all the results obtained, we would like to make some concluding 
remarks. 
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In general, coefficients aj can be computed symbolically by hand and 
only numerically by Mathematica (see Section 3.3); at present we are 
not aware of the possibility to overcome this problem. 

For the known case of spherical barotropic vorticity equation, values 
of coefficients ai coincide with known form the literature for all triads 
but three. These 3 triads, though satisfying resonant conditions, are 
known to be special from the physical point of view in the following 
sense (see [16] for details). Though resonance conditions are fulfilled 
for the waves of these triads, they, so to say, do not have a place in the 
physical space to interact and their influence (if any) on the dynamics 
of the wave system has to be studied separately from all other waves. 
Our results might indicate that also the coefficients of these triads 
have to be defined in some other way compare to other resonant triads. 
For instance, another way of space-averaging has to be chosen. 

The results of Section 3.4.2 show that analytical formulae given in [15] 
for ai are not correct. 

The results of Section 4.3 show a crucial dependence of the number 
of solutions on the form of the boundary conditions. In particular, 
some boundary conditions (for example, {Lx,Ly) = (11,29)) yield no 
solutions which is of most importance for physical applications. From 
the mathematical point of view, an interesting result has been observed: 
in all our computations (i.e. for m,n < 300) indexes corresponding 
to non-empty classes, turned out to be odd. It would be interesting to 
prove this fact analytically because if it keeps true, we can reduce the 
computational time. 

The algorithm presented in Section 4 has been implemented before 
numerically in Visual Basic, and our purpose here was to show that it 
works fast enough also in Mathematica. The algorithms presented in 
Section 3 and Section 5 have never been implemented before, the whole 
work is usually done by hand and some mistakes as in [T5j are almost 
unavoidable: it takes sometimes a few weeks of skillful researchers to 
compute interaction coefficients of dynamical systems for one specific 
wave system. 

All the algorithms presented above can easily be modified for the case 
of a 4-term mesoscopic system. The only problem left is a procedure 
to establish all non-isomorphic topological elements for a quadruple 
graphs, similar to the procedure given in [17] for a triangle graphs. 
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The structure of quadruple graphs is much more comphcated while 
some mechanisms of energy transfer in the spectral space do exist [2S] 
that are absent in 3-term mesoscopic systems. A complete classification 
of quadruple graphs is still an open question but in a given spectral 
domain it can be done directly (a very time consuming operation). 

• We have developed a Web interface for the presented methods, which 
turns the implementations from only locally available software to Web- 
based services that can be accessed from any computer in the Internet 
that is equipped with a Web browser. The presented implementation 
strategy is simple and based on generally available technologies; it can 
be applied as a blueprint for a large variety of mathematical software. 
In particular, the results are not bound to the current Mathematica 
implementation but can be adapted to any other computer algebra 
system (e.g. Maple) or numerical software system (e.g. MATLAB) of 
similar expressiveness. 

• At present, an explicit form of eigen-modes ([5]), ([6]) is used as one of 
the input parameters for our program package. Theoretically, at least 
for some classes of linear partial differential operators and boundary 
conditions, computing eigen-modes can also be performed symbolically 
basing on the results in [26]. If this were done, not an eigen-mode but 
boundary conditions would play role of input parameter. 
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